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This is the first paper in a series aimed to implement boundary conditions consistent with the 
constraints' propagation in 3D unconstrained numerical relativity. Here we consider spherically 
symmetric black hole spacetimes in vacuum or with a minimally coupled scalar field, within the 
Einstein-Christoffel (EC) symmetric hyperbolic formulation of Einstein's equations. By exploiting 
the characteristic propagation of the main variables and constraints, we are able to single out the 
' only free modes at the outer boundary for these problems. In the vacuum case a single free mode 

exists which corresponds to a gauge freedom, while in the matter case an extra mode exists which is 
associated with the scalar field. We make use of the fact that the EC formulation has no superluminal 
^ , characteristic speeds to excise the singularity. 

We present a second-order, finite difference discretization to treat these scenarios, where we 
' implement these constraint-preserving boundary conditions, and are able to evolve the system for 

essentially unlimited times. As a test of the robustness of our approach, we allow large pulses of 
gauge and scalar field to enter the domain through the outer boundary. We reproduce expected 
results, such as trivial (in the physical sense) evolution in the vacuum case (even in gauge-dynamical 
simulations), and the tail decay for the scalar field. 

m 
o 

I. INTRODUCTION AND OVERVIEW 

In the initial value problem of Einstein's equations one decomposes the original system of equations (given by 
G a b = 87rT a 6) into two distinct sets: one set consisting of evolution equations (involving time derivatives of the main 
variables) and the other consisting of constraint equations (which do not involve time derivatives). There are infinite 
possible ways of achieving this decomposition, and the resulting systems are known by a variety of different names 
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I (e.g. ADM, characteristic-Bondi and conformal-Einstein approaches, to name just three). However, depending on the 
^jjy character of the hypersurfaces used to foliate the spacetime, these formulations can be seen as belonging to different 
groups: the group of Cauchy formulations (which require a spacelike foliation), that of characteristic formulations 
(having a null foliation) , or a more generic group (where the foliation's leaves need not have a fixed specific character) 
k> , (see, for instance, 

Irrespective of the group and restricting to the initial value problem (where the problem is boundary- free), the 
state of the system is defined at an initial hypersurface So and the solution to the future of So is obtained through 
the evolution equations. However, if the system of equations involves constraints, one might opt to employ only a 
subset of the evolution equations supplemented with enough constraint equations to solve for the field variables. This 
strategy is referred to as constrained evolution, as opposed to free evolution (where only the evolution equations are 
employed). In the particular case of Einstein's equations, a straightforward manipulation of the Bianchi identities 
demonstrates that either strategy produces, at the analytical level, the same solution. Therefore, one need not deal 
with constrained evolution (which usually requires solving elliptic equations) and the more direct approach of free 
evolution can be safely employed. In the numerical realm, however, the picture is more complicated (for definiteness, 
from now on we concentrate on the case of Cauchy evolution). On one hand, employing constrained evolutions might 
represent a significant computational overhead as it usually involves solving elliptic equations at each time step; for 
this reason constrained evolutions have been, for the most part, avoided beyond the two-dimensional case. On the 
other hand, free evolution in numerical implementations (which only evaluate the constraint equations to monitor 
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the quality of the implementation) display violation of the constraints. At early times, these violations are consistent 
with the truncation error Q , but as time progresses the observed violations often grow quite rapidly. 

The picture gets even more complicated in the presence of boundaries. The violation of the constraints in uncon- 
strained numerical evolutions frequently grows without bound. Possible reasons, besides the ones already present 
in boundary-free models, are constraint-violating modes introduced by standard boundary conditions (which might 
drive instabilities, and in any case do introduce spurious errors that should be avoided). If the case under study 
corresponds to cosmological scenarios or short term evolutions, the aforementioned problem should not be worrisome. 
In the cosmological case one has no boundaries (or rather, periodic boundary conditions are specified, which trivialize 
the specification issue) , for solutions which are restricted to the future domain of dependence of the initial data no 
boundary conditions need be specified, and for short-term evolutions the errors introduced by approximate boundary 
conditions are restricted to a small region. There are, however, many non-cosmological important problems that 
require long term evolutions and hence the need to prescribe consistent boundary conditions; among these, collapse 
scenarios and black hole spacetimes. 

In the one-dimensional case, approximate boundary conditions that minimize the influence of errors introduced at 
the boundary have been presented ||. This is achieved by placing the boundaries far away and treating a region 
close to them in a special way (to control spurious reflections). Unfortunately, these techniques are not as effective 
in higher dimensional scenarios, as the computational cost involved in placing boundaries far away is excessive. 
Furthermore, even when there were enough computational resources, it would be preferable to make use of them 
to achieve finer resolved simulations, rather than to push boundaries farther away. It is therefore of considerable 
importance to formulate relatively inexpensive boundary conditions whose associated errors do not depend on the 
boundaries' locations, minimize spurious reflections, and guarantee that constraint-violating modes are not introduced. 

Essentially, one aims to provide boundary conditions such that there is not only a unique solution to the evolution 
equations, but there is also preservation of the constraints throughout the computational domain. In a strongly 
hyperbolic formulation of gravity (see for reviews on hyperbolic techniques applied to Einstein's equations) the 
first requirement amounts to giving boundary conditions only to the characteristic modes that enter the computational 
domain at any given boundary. Satisfying the second requirement is considerably more involved and has only very 
recently started receiving attention. In the analytical realm, well poscdncss of the initial boundary value problem 
(for a particular system) has been established in g , which sheds light on the physical understanding of the issue and 
shows that, at least in a by-case basis, the problem might be analytically tractable. In the numerical realm, several 
efforts (restricted to different scenarios) have illustrated the advantages of providing boundary conditions through 
the use of constraints [||. It is precisely this problem that we want to address here within the context of Cauchy, 
unconstrained evolution. 

The present work aims to contribute to the area with the particular motivation of addressing issues proper of 
numerical implementations. Although in this first paper we restrict our analysis to the spherically symmetric case 
and a particular way of reformulating Einstein's equations, the procedure is straightforward to implement within any 
hyperbolic formulation of gravity and in three dimensions (3D). 

We present a detailed discussion of the analytic treatment for the boundary conditions and illustrate its benefits 
with a simple numerical implementation which accurately evolves a spherically symmetric spacetime, stationary or 
dynamical (the dynamics either corresponding to gauge modes in the vacuum case or to true dynamics in the case 
of a minimally coupled scalar field), for unlimited periods of time and without sign of instabilities. Throughout this 
work we minimize the use of techniques specifically developed for treating spherically symmetric spacetimes, both 
analytical and numerical, to expedite the generalization to the 3D case. 

Our starting point is the assumption that we are dealing with a first-order, quasilinear strongly hyperbolic formula- 
tion for Einstein's equations |]. In the spherically symmetry case such a system can be written asii = Au' +l.o., where 
u is an array of variables, dot and prime indicate time and spatial derivatives, respectively, A (called the principal 
pari) is a diagonalizable matrix that might depend on u and the spacetime coordinates, but not on derivatives of u, 
and l.o. stands for lower order terms, i.e. terms that do not have derivatives (of any kind). So, we have a system of 
equations for the main variables^, 



1 In 3D, a generic first-order system can be written as ii — y\ A^diU + l.o.; defining the symbol P := A l uii, such a 
system is said to be symmetric (or symmetrizable, the terminology in the literature is not uniform) hyperbolic if P| w =i| can be 
symmetrized, by a smooth transformation independent of Co := (u\,uJ2, W3), while it is said to be strongly hyperbolic if it can 
be symmetrized (the symmetrizer is not required to be independent of Co), and its eigenvalues are real. In spherical symmetry, 
a strongly hyperbolic system is symmetrizable (the converse statement is always true), but this is an artifact of one dimension 
(ID). 

2 The word main will be used, when there is possibility of confusion, to differentiate between the evolution equations for say, 
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u = Av! + l.o. (System I) , 



(1) 



and evolution for the constraints, which we assume are also strongly hyperbolic, 



u c = A c u' c + l.o. (System II) . 



If one were interested only in system I, one would give initial data for the variables that form u, and boundary 
data only to the ingoing characteristic modes^J. These characteristic modes (eigenvectors of A) are, depending on 
which boundary one is dealing with, the ones that are travelling to the left (positive eigenvalues) or right (negative 
eigenvalues) . 

A unique solution to system II is fixed, similarly, by giving initial data to u c and boundary conditions to the ingoing 
modes of the system. Since this is supposed to be an homogeneous system (which is the case in Einstein's equations 
) , the identically zero solution is obtained by providing zero as initial data (u c = 0) and zero boundary conditions 
to the eigenmodes of A c entering the domain. The crucial point here is that, in unconstrained evolution, initial data 
and boundary conditions for u c follow from that of u. Thus, one has to provide initial and boundary data for u 
such that they imply zero initial data and boundary conditions for u c (i.e. the data are consistent with all Einstein's 
equations). The first step, namely, providing initial data that satisfy the constraints, is customarily well satisfied as 
the constraint equations are employed in the determination of data for u (see [^JQ] and references cited therein). It is 
the second step that is usually neglected. Here we concentrate on this problem with the goal of providing consistent 
data, whose associated error neither depends on the location of the boundaries nor on the total evolution length; 
rather its error will agree with the overall truncation error of the global implementation. 

The main difficulty with this second step is that the constraints (and, therefore the eigenmodes of system II as well) 
involve spatial derivatives of the main variables. On the other hand, by controlling the (ingoing) characteristic modes 
of system I at the boundaries, one does not control the spatial derivatives needed to set the ingoing constraint modes 
to zero. One way around this (and the one we use here), is to trade spatial derivatives for time derivatives using 
system I. Consequently, enforcing the constraints at the boundary amounts to controlling some time derivatives. We 
will make this construction explicit later but, before proceeding further, some comments are appropriate (both issues 
certainly deserve further analysis): 

• To our knowledge, there is no rigorous result guaranteeing that in any strongly hyperbolic formulation of 
Einstein's equations this trading of spatial for time derivatives can be done. 

• When performing the trading, one ends up with some conditions at the boundaries on some time derivatives of 
the main variables. One must find out whether these conditions can be fulfilled by controlling only the ingoing 
modes of system I (at the numerical level, conditions on the time derivatives are enough for the application of 
the method of lines |^]). Again, we know of no proof showing that this should be possible in a general case. 

We here show how this inversion can be performed in the ID case, and leave for a future paper a similar analysis in 
3D linear gravity p0| . 

When dealing with black hole spacetimes, boundary conditions customarily refer to the outer boundary conditions, 
but one might also have inner boundaries. These constitute "holes" or "excised" regions from a given computational 
domain; the most widely considered are those where the excised region has been chosen so as to remove the singularities 
when dealing with black hole spacetimes (assuming the validity of cosmic censorship). In this case, a region inside 
the black hole is excised from the computational domain (Unruh, cited in ]TT[|). This excision strategy introduces an 
inner boundary which, if chosen inside the black hole, should leave the region outside the event horizon unaffected. 
Mathematically, the realization of this idea is ensured by employing a hyperbolic formulation with no superluminal 
characteristic speeds and, in particular, where all eigenvalues describe modes propagating towards the inner boundary^. 
The non-triviality of excision if one does not have a formulation with these properties is discussed in one of the 
appendixes. 

In recent years, a number of re-formulations of Einstein's equations with physical characteristic speeds have been 
presented [O-pTj]: Here we make use of one of them in our analysis, the so-called Einstein-Christoffel (EC) system ll^] 




(note that any other would have been equally well suited - our choice is simply motivated by a comparison with available 
results (15)). 



the three-metric and extrinsic curvature (and possibly extra variables) and evolution for the constraints. 
3 By ingoing we refer to those modes entering the computational domain with respect to a given boundary. 
4 For a discussion of different approaches towards application of excision techniques see llf and references cited therein. 
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The organization of this paper is the following: in section II we present the equations for the EC system describing 
a spherically symmetric massless scalar held minimally coupled to gravity. In section III we describe how to give 
boundary conditions that preserve the constraints, while in section IV we present our numerical method to test the 
procedure. Numerical results are included in section V, and are divided into three classes: evolution of stationary 
slicings of a Schwarzschild black hole, evolution of dynamical slicings of the same spacetime, and evolution of a massless 
scalar held interacting with a black hole. In all cases the simulations can be followed for unlimited times, even using 
very moderate resolutions. We demonstrate the accuracy of our method by monitoring the mass of the black hole 
in the vacuum case and reproducing the known tail decay in the scalar held case. A particularly strong test of the 
robustness of the approach is presented by letting strong gauge or scalar held pulses enter the computational domain 
through the outer boundary (compare with standard outer boundary treatments, where conditions are obtained by 
requiring the geometry at the outer boundary to be close to flat or Schwarzschild spacetimes, i.e. quite the opposite 
from what we do here). 

II. MAIN EVOLUTION EQUATIONS AND PROPAGATION OF THE CONSTRAINTS 

As often is the case in hyperbolic formulations of Einstein's equations, the EC formulation uses "exact" (i.e. arbitrary 
but a priori specified) shift and exact densitized lapse (defined by a = Ng~ x / 2 , where g is the determinant of the 
three-metric and N is the lapse). In the present work, and for the sake of maintaining a uniform notation, we will 
follow as much as possible the conventions of |lE|| . For example, we write a :— ar 2 sin 9, and 

ds 2 = -N 2 dt 2 + g rr (dr + fldt) 2 + r 2 g T {d6 2 + sin 2 8d(j) 2 ) , 
K i6 = K rr dr 2 + r 2 K T (d0 2 + sin 2 8d<j) 2 ) , 

where all fields depend only on (t,r). Since the EC formulation is a first-order reformulation of Einstein's equations, 
besides the three metric and extrinsic curvature, further variables (which basically contain information of spatial 
derivatives of the three-metric) are needed. In our present case of spherical symmetry, only two new variables require 
introduction, f rrr and f r x, defined by 

f - 9 rr i 4 9rrfrT 
Jrrr — _ T j 
2 9T 

f 9t , 9t 

JrT = -7T H • 

2 r 

Additionally, we plan to study a scalar field ^ minimally coupled to the geometry. In order to re-express the equation 
that governs this field, 

g ab V a V b y = , 

as a first-order hyperbolic system, we introduce two new variables: 

n := -(/?*' - *) , 

a 

$:=*'. 

The evolution equation for W decouples from the rest, in the sense that one has a closed system for the set of eight 
variables (g rr , gr, K rr , Kt, frrr, /rT,n, $) which one can solve for, and afterwards obtain 'P. Because of this, and 
following common practice, we will drop \? from the system. 

Thus, the main evolution equations, up to principal part (the complete expressions are listed in the appendix), are 



9rr = Pg'rr + 1-0. , (2) 

9T = Pg' T + l.o. , (3) 
N 

K rr = PKr frrr + l-O. , (4) 

9rr 

N 

K T = pK' T f rT + l.o., , (5) 

9rr 

frrr = (3f r „ ~ N K' rr + l.O. , (6) 
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frT = Pf'rT ~ ^ K' T + 1.0. , (7) 

N 

n = pn! + i.o. , (8) 

Qrr 

$ = - Nil' + l.o. . (9) 

The characteristic modes and eigenvalues determined by the system play a crucial role in our boundary treatment. 
These are (note that the modes with speed j3 propagate along the timelike normal to the foliation, while the other 
modes propagate along the light cone), 

ui = g rr (v 1 = (3) , (10) 
u 2 = g T (v 2 = (3) , (11) 

U 3 = K rr - grr /2 frrr ("3 = P + &9t) , (12) 

u A = K T - grr 1/2 frT (vi = (3 + ag T ) , (13) 

U 5 = K Tr + 9rr /2 frrr («5 = P ~ «5t) , (14) 

u 6 = K T + grr /2 frT (v 6 = (3 - ag T ) , (15) 
u 7 = n + 5 - 1 / 2 $ (v 7 = p-ag T ), (16) 
u 8 =H-g- 1 / 2 ^ (v 8 =(3 + &g T ). (17) 

If we were dealing with a constraint-free system described by equations (^[^), a unique solution would be fixed by 
providing boundary conditions for the characteristic modes that are incoming at the boundaries and data on the 
initial hypersurface. However, we know that a solution to equations (^j^) is a solution of Einstein's equations if and 
only if the following four constraints are also satisfied: 

C:= J^_ 1 + W2 + 7^_A^ _W^ + KT\ + ^tf 0> (18) 
g rr gr ^ z g T g rr gT \r 2g T g rr J g T \ g rr 2g T J ig rr 4 

C r := — + - — (— + — +$n = 0, (19) 



9T rg T gr \ g rr gr 

C rrr . — g rr -\- — 1f rrr — 0, (^0) 

g T 

C rT — - VrT = . (21) 

r 

The first two are basically the Hamiltonian and momentum constraints, respectively, while the other two correspond to 
the definitions of the extra variables that make the system first order with respect to spatial derivatives. As mentioned, 
it is common practice to choose consistent initial data for Einstein's equations by solving these constraints; however 
the constraints have been examined in a limited number of cases to provide consistent boundary data. The main 
purpose of this work is to provide further indications that constraints should be looked at more closely when dealing 
with boundary conditions. First, note that these constraints are defined in terms of the main variables, and a solution 
of the main evolution equations completely determines them as functions of spacetime. In particular, one can obtain 
the time evolution of the constraints by: (i) taking time derivatives of the right hand side (r.h.s.) of equations (p^-pl]); 
(ii) replacing the time derivatives by the r.h.s. of equations (|[|§); and, (iii) re-expressing the main variables in terms 
of the constraints and their spatial derivatives. In the present case, these are 

N 

C = (3C C' r + l.o. , 

Qrr 

C r = -NC + (3C' r + l.o. , 
C rrr — (3C, rrr -f- l.o. , 
C rT = (3C' rT +l.o.. 

We can also calculate the characteristic modes and eigenvalues of this system, obtaining: 

Ci =C + g- 1 ' 2 C r {vt = 0-&ffr), (22) 
C 2 = C-g-V 2 C r (v c 2 =[3 + a g T ) : (23) 

C 3 = Crrr («§ = P), (24) 

C 4 = CVt « = /3)- (25) 
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In the next section we exploit this information about the characteristic structure of our system in order to set up 
boundary conditions that will preserve the constraints at the boundaries. 



III. CONSTRAINT-PRESERVING BOUNDARY CONDITIONS 



In any system where the characteristic modes and eigenvalues are known, the data required at a given boundary 
are straightforwardly assessed by examining the eigenvalue of each characteristic mode. In our present case, these 
eigenvalues are Ai = ft — agr, A2 = /3 + olQt-, and A3 = f3. In the case where the shift is exact, one can a priori 
guarantee which sign A3 will have. Throughout this paper we shall take as positive, since in that way we will be 
able to reproduce stationary known slicings of the Schwarzschild space-time^], and to evolve dynamical spacetimes 
as well. The positivity of the shift implies that A2 is also positive. On the other hand, the sign of Ai depends on 
the solution and, thus, cannot be controlled a priori. However, for typical stationary slicings of Schwarzschild it is 
negative (positive) outside (inside) the black hole, and zero at the horizon. By continuity, the same will hold for 
(perhaps slight) distortions of these slicings. In our simulations we check numerically that this is indeed the case. As 
we shall see, even on highly distorted spacetimes this condition remains satisfied. Figure 1 shows a schematic diagram 
for the characteristic modes of a Schwarzschild black hole. 



FIG. 1. Schematic diagram for the characteristic speeds of the different modes in a stationary slicing of Schwarzschild black 



Our inner boundary is always inside the black hole, and no conditions are needed there. The usual procedure to 

make sure that the inner boundary is inside the hole is to track the apparent horizon. In spherical symmetry this is 

1/2 

particularly simple since its location is defined in terms of algebraic combinations of the main variables (f r T—grr Kt = 
0). Our actual procedure is to choose the initial data such that the inner boundary is inside the apparent horizon, 
and then numerically monitor its position during evolution. In our simulations the horizon never touches the inner 
boundary, so we do not need to move it, i.e. its grid location is fixed. 

From the discussion in section II, we know that at the outer boundary (which in our case also has fixed grid location) 
we need to specify u\, U2, U3, 114, from the characteristic structure of the main equations. However, since the system 
is constrained, we are not free to choose these arbitrarily (which would be the case in an unconstrained system), 
but rather we must do so such that C2 = C3 = C4 = are satisfied (Note that C\ is outgoing, therefore nothing 
special needs be done for it, nor should it be, otherwise we would overdetermine the constraint system of evolution 
equations.). This procedure, coupled with initial data satisfying the constraints, and the fact that the constraints are 
propagated with a quasilinear homogeneous system, will ensure that they are preserved everywhere. 

We now make explicit how this procedure is implemented in our present case. We start by discussing how to enforce 
C3 = C4 = at the (outer, from now on) boundary. Writing down the constraints explicitly in equations ( p4|p5| ), 



5 Kerr-Schild, Painleve-Gullstrand, full harmonic and time harmonic slicings of Schwarzschild have positive shifts, see Jl4| for 
the explicit form of these metrics in the EC formulation. 




hole. 



G 



these conditions are 



g'rr + 



JrrfrT 
9T 



9t 



2g T 



2,f rrr 0, 

- 2frT = . 



Using equations (||||), these can be rewritten as 



g rr = 2g rr /3' — 2NK rr — 
g T = -2NK T + 20f rT . 



, 9rrfrT 



(3 + 2f rrr (3 



(26) 
(27) 

(28) 
(29) 



Note that since g rr and gx correspond to ingoing modes, we are free to choose them such that they satisfy Eqs. ( p8| , f29| ). 
As was mentioned, one does not need to solve these two differential equations at the boundary, since time derivatives 
are all that is required for the time update in the method of lines, which we use for the time update (see section IV). 
Thus, we simply use the expressions (28 2^) as the right hand side of the corresponding field variable at the boundary 
points (and similarly for f r j- and Kj-, see Eqs. fl3l| , |32| ) below). 

Lastly, C2 = needs to be enforced at the boundary. Using its definition, we have 

-JrT - K T + I* = • 



V9r 

Trading spatial for time derivatives, this in turn gives 

1 



>9r 



-J, 



rT 



K T + l.o. = 



(30) 



Next observe that, by definition, Kt = Kt{ua,uq, 9rr) and f r T — frT{uA,u§, 9rr)- Since g rr at the boundary is 
readily known from equation J2S|) , Eq. ( |30| ) fixes the ingoing mode 114 (up to now free). From this, the definition of 
U4 and us, and Eq. (Eq), we end with (recall that uq is outgoing, so it does not need boundary conditions) 



1/2 

i " rr ■ 1 ; 

JrT = ~~^~ u f> + Lo - 

K T = 7^U 6 + l.O. 



Similarly, the ingoing mode U3 is completely arbitrary. From it and the outgoing mode u§ we have 



K rr = ~(u 3 + u 5 ) , 
1/2 

frrr ^ (^5 



Finally, the ingoing mode ug, is also arbitrary. From it and the outgoing mode u-j we have 



n = 2^ U7 + U8 "> ' 

1/2 

■i 9rr / x 
<P = —^-{U7 - U 8 ) ■ 



(31) 
(32) 

(33) 
(34) 

(35) 
(36) 



Equations (|2^,^,^T|,^2 33 3^ , |35| , |36|) completely determine the boundary conditions for our eight variables; U3 and 
us are the only free modes. In the vacuum case, u$ = and U3 describes a gauge mode. 

In the next section we discuss our numerical implementation and results. One comment is in order before that: 
there is an additional constraint, C ma tter = $ — but by repeating the treatment described above one can see that 
this constraint fixes the boundary condition for 'J: Since we are not evolving 5', we do not need to care about this. 
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IV. NUMERICAL SIMULATIONS 



To implement the proposed boundary treatment strategy, we have chosen a straightforward second-order dissipative 



method of lines. In particular, it was not necessary to employ techniques such as causal differencing |16|, upwind 
discretization or any special way of treating the Lie-derived terms in Eqs. (|^-^|) (however, these can be readily 
incorporated as well). Spatial derivatives are discretized with second-order centered differences plus fourth-order 
dissipation as discussed in Jl7| , while for the time integrator we use second-order Runge-Kutta (the dissipation added 
is, indeed, needed; otherwise, this method is unstable even for a simple scalar wave equation; see, e.g., ||). 

Our uniform grid structure consists of points i = ... N, with grid spacing Ar = L/N, where L — r outer — 
and we implement derivatives according to standard formulae 

Af -> AD f - ^{ArfD+D+D.D^f , 

/ j-, r\ fi+l ~ fi—1 



' inner j 



(D+fh 
(D-f)i 



2Ar 

fi+l ~~ fi 

Ar 

fi — fi-1 

Ar 



To evaluate derivatives at the boundaries, one either resorts to one-sided derivatives (which require different algo- 
rithms applied at boundary points) or introduces ghost zones, which are artificial points beyond the boundaries where 
field values are defined via extrapolation (one can choose the extrapolation order such that the answers from the 
different approaches are exactly the same) . For convenience we chose the latter approach with only one ghost-zone for 
each boundary. Field values at these ghost zones are defined via third-order extrapolations and the same derivative 
operator is applied Vi, i = 0..N. 

Since our inner boundary is always inside the black hole, we extrapolate all variables at the ghost zone point 
(* = -!)• 

At the outer boundary ghost zone (i = N + l), the outgoing modes U5 and U7 are also found by extrapolation, while 
the ingoing modes U3 and us are set as arbitrary functions of time. Next, equations ( |2§| , |29| , |3l| , |32| ) are integrated, at 



each time step, with second-order Runge Kutta. This, plus equations ([53|,[34 35 3^), completes the treatment for the 
outer boundary. 

As a side note, it is worth mentioning that since we use only one ghost zone for each boundary, fourth-order 
derivatives cannot be obtained at the points i = 0, N; thus, no dissipation is added to these points (we set o- i=0 = 

0- i=N = 0). 

In all the cases discussed in this paper we have checked second-order self-convergence for the eight main variables, for 
the constraints and mass, as well as convergence with respect to the analytic solution, whenever this is available. The 
Courant-Friedrich-Levy factor A = At/ Ar is set to 0.25, and the dissipation factor is ai — 0.6/16 (for i — 1..N — 1). 

A. Vacuum evolutions 

1. Stationary slicing of a black hole 

The simplest test of a black hole spacetime consists of reproducing known stationary slicings of a Schwarzschild 
black hole. For this test we give the corresponding known values both to the initial data and to the 113 mode (its=0 for 
vacuum) at the outer boundary. We concentrate here on Painlcvc-Gullstrand slicings, but we have obtained similar 
results using Kerr-Schild slicings. The code runs for unlimited times, even with resolutions as coarse as Ar = M/6 
(with M the mass of the black hole) . Figure 2 displays the L2 norm of the Hamiltonian constraint for three different 
resolutions and L = 9M (r 'inner — M), the L2 norm of a grid function / being defined as 



1 7V ~ 1 \ 

i=0 / 



1/2 



l/h 

There is no growth in the Hamiltonian constraint, and the same holds true for the other three constraints. 
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Ar-M/16 

■-■ Ar-M/32 




0.0 001 1 1 ' — 1 — ' 1 1 — 1 < — < — 1 

10 100 1000 10000 

t[M] 

FIG. 2. 1/2 norm of the Hamiltonian constraint, for the evolution of a Painleve-Gullstrand black hole. In these runs the 
outer boundary is at r = 10M. 
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FIG. 3. Z/2 norm of the relative mass error \8M\ = \M — M\/M, for the same runs shown in figure 2. After one thousand 
crossing times, the relative error for the coarsest resolution is around ten percent. 



Figure 3 shows, for the same run, the L<i norm of the relative error, |5M|, in the mass function M, denned as 

\5M\ :=L\ M -M\ , 

where M, the Misner-Sharp mass ]l8| , is a gauge invariant |]l9| definition of the ADM mass M in spherical symmetric 
vacuum: 

+ lL(K^-^S\ . (37) 

9T \ 9rr J _ 

That is, each of the terms in the r.h.s. of Eq. ( |37| ) might be a function of t and r, but Einstein's vacuum equations 
in spherical symmetry imply that, at the continuum, the r.h.s. is a constant (equal to M), both as a function of space 
and time. When we compute SM we evaluate M. using the numerical values of the right hand side of Eq. ([57|), while 
for M we use its analytical value. 

From figure 3 we see that, as opposed to the constraints, there is a drift in the mass (this effect is expected, as 
second-order errors accumulate after many iterations) which, after a couple of hundred masses, is essentially linear in 
time. In figure 4 we plot M. /M as a function of radius, for different times and Ar = M/8; the error is spread all over 
the domain and the mass decreases monotonically in time, which suggests that the black hole is moving through the 
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FIG. 4. Numerical value of the Misner- Sharp mass, scaled by its analytical value, as a function of radius, for different times. 
The monotonic behavior is explained by the fact that the black hole "moves" through the grid. 



grid (in the sense that r — > r + t x v num ). We have verified that this is indeed the case (but, still, v num = O (Ar 2 )) 
by computing the (grid) location of the apparent horizon. For example, with Ar = AT/8, after 10,000M the horizon 
has moved from r^ = 2M to = 1.75M, which accounts for the roughly ten percent error that the mass has at that 
time in figure 3: 

SM w l-r h /(2M) re 0.1 

We have tried different positions for both boundaries, with the inner one always inside the black hole and the outer 
one ranging from 10M to 1, 000 M , and neither the stability nor the convergence of the simulations depended on such 
positions. In figures 5 and 6 we show the analogue of figures 2 and 3; with the outer boundary placed at 50 M , while 
still running the code up to 10 3 crossing times (i.e. t = 50, 000M). 
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FIG. 5. Same as for figure 2, but now with the outer boundary at r = 50M. 
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FIG. 6. Same as for figure 3, but now with the outer boundary at r = 50M, as in the previous plot. 



2. Gauge pulse falling into black hole 

As discussed in Section III, the mode 113 at the outer boundary is completely arbitrary. In this next test we let a 
gauge pulse enter the domain through the boundary. We thus use initial data corresponding to a Painleve-Gullstrand 
slicing of a black hole of mass M, and fix M3 at r outer by superposing, on top of the Painleve-Gullstrand value, a 
Gaussian pulse described by 



,PG 



1 + Ae -(*-*o) 2 



(38) 
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FIG. 7. The radial component of the metric, for a gauge pulse entering the domain through the outer boundary. The 
amplitude's pulse grows as it approaches the inner boundary, falls into the black hole, and finally the metric settles down to a 
stationary one. 



In hgure 7 we present snapshots of the g rr component of the metric as a function of radius for different times. The 
width of the pulse is a — 2M, and it is centered at to = 5M. The amplitude is quite large, A — 1, corresponding 
to a 100% "disturbance" of the stationary value. The resolution employed is Ar = M/8, and the domain extends 
from ri nner = M to r ou t er = 30M. The pulse grows as it approaches the inner boundary (which is not surprising, 
since it happens even for a linear scalar equation propagating in a Schwarzschild background due to the effective 
non-trivial potential) , g rr having a disturbance of more than 500% that of the stationary value when it "crosses" the 
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inner boundary. After the pulse falls into the black hole, the metric components gradually settle down to stationary 
values and the code runs for unlimited times. 

Since we are dealing with a vacuum spherically symmetric spacetime, the resulting spacetime must be a (dynamical) 
slicing of Schwarzschild. One way of corroborating this is to follow the Misner-Sharp mass. As can be seen in figure 
8, where \5M\ is shown as a function of time for different resolutions, the numerical value of the mass does converge 
to this analytical prediction. Just as in the stationary case, the code runs for as long as wanted, and has a similar 
linear drift in the mass at late times. The growth in the errors around t = 30A/ seen in figure 8 are not due to the 
pulse entering the domain through the boundary (this happens at t = 5M) but due to the pulse reaching the inner 
boundary, where all the gradients are steeper, and the huge growth in g rr is rather poorly resolved. Similar growths 
appear in the constraints (see figure 9), but they are still second-order convergent. It is indicative of the power of 
consistent boundary conditions that we can have such a big pulse entering the domain through the boundary, causing 
the spacetime to be so dynamical, that the code is stable and that the mass error is small (notice that even for a 
resolution as coarse as Ar = M/8 this error, while the pulse travels trough the domain, does not exceed one percent). 




t[M] 

FIG. 8. Mass error for a gauge pulse "falling" into the black hole. The peaks around t — 30M are caused by the pulse 
reaching the inner boundary. After the pulse falls into the black hole the spacetime settles to a stationary one, and at late 
times the mass has a linear drift, just as in the stationary evolutions. 
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FIG. 9. Z/2 norm of the Hamiltonian constraint for the same run shown in figure 8. 

In equation ( |§8| ) we wrote the boundary condition in that way (i.e. the Painleve-Gullstrand value times a time 
dependent function) in order to have a measure of how large the pulse is compared to its stationary value. In the 
next simulations we will have a more physical idea of how strong these pulses can be in our code. 
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B. Scalar field coupled to the black hole 



As a last test of the approach we consider the case of a minimally coupled scalar field which also enters the 
computational domain through the outer boundary. As mentioned, only one scalar field mode is freely specifiable at 
the outer boundary. It is this mode which is chosen to describe an incoming pulse with amplitude A and compact 
support in time as 

u 7 (t) = A(t - t/) 4 (i - t F f sin(Trt) , 

if t G [t],tp], and Uf(t) = otherwise. We choose U3 = 0, and for the initial data we take that of a Painleve-Gullstrand 
black hole of mass M. Note that one advantage of this is that one can perform non trivial simulations without solving 
the constraints initially, since one can provide, as we have done, a known solution of the constraints as initial data, 
and introduce the non trivial dynamics through the boundary (where all the treatment is algebraic), see also pofl . 

We have tried different amplitudes and time range for the pulse, obtaining stable discretizations in all cases. In 
order to illustrate the robustness of the approach we here show two examples. In the first one we show how the 
method can cope with considerably strong pulses (measured by the amount of energy that the black hole accretes), 
and in the second one the tails of the scalar field are computed and shown to agree with the expected result. 



100 

t[M] 



FIG. 10. Change in the area of the apparent horizon, due to a quite narrow (At = 10M) pulse of scalar field injected at the 
outer boundary, with a not very demanding resolution (Ar = M/8). The final mass of the black hole is more than twice the 
original one. 



For the first case we chose A = 8, tj = 0, ip = 10M, Ar = M/10 and the outer boundary at r = 50M. Figure 
10 illustrates the area of the apparent horizon as a function of time. Initially the area is M and, as time progresses, 
it increases by 270%, until it reaches a stationary state describing, as expected, a Schwarzschild black hole of larger 
mass. 

In the second case, in order to resolve the tails accurately, considerably finer resolution is needed. For this reason, 
we chose Ar = M/50 and place the outer boundary at r = 100M. In order to measure the tails, we place four 
different observers: at r = 50, 70, 90 and 100M. The decay rate (of the form $ cx t~ n ) found by these observers 
is n = —2.89, —2.97, —2.94 and —3.04, respectively. These are in excellent agreement with the expected value of 
n = —3 pl| ]. Furthermore, the clean treatment of the outer boundary allows for accurate measurements even at 
the last point of the computational domain. Past works, which resorted to approximate boundary conditions, have 
observed that the boundary influenced the results when the observers were placed close to it. 

Finally, we should mention that we have tried with different kinds of time dependent boundary conditions for the 
modes U3 and uj, and the results presented here do not depend on the particular choice made. Not only can we define 
pulses that last for a finite amount of time, as we have shown here, but we can define a time dependency that lasts 
forever (say, U3 = Asnx(u)t)). In this sense, the large growth in the mass of the black hole mentioned above is what 
one can do with a pulse of width 10Af. In fact, one can increase arbitrarily the mass by using wider pulses, i.e. by 
slowly injecting scalar field energy through the boundary. 
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FIG. 11. Tail decay for the (time derivative of the) scalar field. The observer is located at the last grid point (r = 100M) 
and measures a tail decay of <E> oc t -3 04 , in good agreement with the expected decay of t~ 3 . This is so mainly because of the 
clean treatment for the boundary. 



V. CONCLUSIONS 



In the present work we have presented additional support for the need to consider the constraint equations when 
defining boundary conditions. Although approximate boundary conditions (which are simpler to implement) can 
certainly be arranged to produce reliable simulations, for long term evolutions they have some important disadvantages: 

• The inherent error introduced might accumulate, spoiling the physical outcome: For instance, even in the 
particular case of a binary black hole collision, the radiation is expected to be around 5% of the total mass of 
the system. Therefore, even small errors which accumulate over time can account for a significant percentage of 
the radiated energy. Since the outcome of some of these simulations will be used as templates for gravitational 
wave data analysis p2| ], one should minimize all possible waveform contaminating sources and diminish both 
amplitude and phase errors. 

• Although the inherent error of approximate boundary conditions can be reduced by placing the outer boundary 
farther away from the sources, this amounts to having to compute the solution on a larger computational domain 
which naturally increases the cost of the simulation. 

• Standard approximate boundary approaches pay little or no attention to the satisfaction of the constraints at the 
boundary. It is hard to conceive that they do not introduce constraint violating modes into the solution. Unless 
the implemented scheme is capable of damping these modes away, they will be present in the computational 
domain. A possible source of instabilities (though not necessarily the only one) in present implementations 
are, precisely, constraint violating modes, hence, it is important to eliminate all possible sources for these 
modes. It might happen that, in certain cases, the observed constraint violations are due to those particular 
implementations. Nevertheless, as the goal is to solve Einstein equations, sources of constraint violations must 
be removed. As mentioned, standard boundary treatments are likely to introduce these violations. 

Furthermore, carefully designed boundary conditions are important in the case of artificial boundaries. These are 
boundaries which occur inside the computational domain where the computational mesh has regions with different 
resolutions (for instance in the case of adaptive mesh refinement) or has been subdivided into sub-domains which are 
treated independently (e.g. when using domain decomposition techniques). As these boundaries are purely artificial 
and introduced for convenience, it is imperative to count with a clean boundary treatment which eliminates (or at 
least minimizes) any spurious reflections or numerical noise. 

Consequently, the search for accurate boundary conditions is of importance in present applications. The results 
presented in this work attest to that effect and suggest a practical and simple way to obtain such boundary conditions, 
which ensures constraint violating modes are absent by the very definition of the approach. Current work is devoted 
to applying the same techniques to 3D problems. 
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APPENDIX A: EVOLUTION EQUATIONS 

The main variables propagate according to 
9rr = Pg' rr + 2g rr (3' - 2ag x J 2 g T K rr 



Wg 



T 



g T = (3g' T - 2ag 1 r i 2 g T K T + 

r 

K rr = f3K' rr - ag- r 1/2 g T frrr ~ a gH?g T - Qg^dl^^frT + ^grr' 1 g\ 12 'a' - 6g T r~ 2 g x J 2 a + 2K rr (3' - 
2g T g 1 r / r 2 a<$> 2 - g T g~ 1/2 aK 2 r + 2g\ l r 2 uK rr K T - 8g- 1/2 af rT frrr + 2g T g„ 3/2 &f 2 rr + 

2gTr~ 1 g~ r 1/2 &f rrr - g T grr /2& ' frrr 

K T = (3K' T - ag T g~}l 2 f rT + 2{3r~ x K T + g T r- 2 g}J 2 a + ag T K T K rr g-^ 2 - gTfrTa'g- 1/2 - 2af 2 T g; r 1 ' 2 
fwr = Pfrrr ~ ag\l 2 g T K' rr - Ag% 2 a'K T + 12g T 1 g^ 2 aK T f rT - kg\ l2 aK T f rrr - g T grr /2 aK rr f rrr - 

10gU 2 aK rr f rT + 3f rrr (3' + g rr p" - olg\l 2 g T K rr + 2r- 1 g T g 1 r l 2 &K rr + Sr^g^aKr + 4a^/ 2 5T $n 

frT = PfrT ~ *9tP '9tK' t + f3' f rT - o! g\' 2 9 T K T + 2g\l 2 &K T frT ~ a 9 ^' 2 K T f rrr g T + 2^0/^ 

$ = /3*' - agU 2 g T n' - g-^&ffrllfrrr + 2agU 2 TifrT + 2r" x ag]' 2 g T U - a'g^grU + 

n = - g- 1/2 &g T & + g- r 1/2 ag T TlK rr + 2ag 1 r /2 nK T - 4g„ 1/2 a$f rT + 2r- 1 g„ 1/2 a gT <S> - g^ ,2 g T <5>a 

APPENDIX B: EXCISION OF THE SINGULARITY 

The issue of excision is a non trivial one, even if dealing with a strongly hyperbolic formulation of Einstein's 
equations, and even ignoring the constraints. The point is that, depending on the formulation, modes can leave the 
black hole. Even though these modes cannot be physical, they would need boundary conditions in any case in order 
to fix the solution. If one, for example, extrapolated all variables in such a situation, one would be implicitly be giving 
boundary conditions that depend on the discretization and that might not have a consistent limit as the grid spacing 
is decreased. Here we show a simple example to illustrate this point. 

We start with a formulation widely used in numerical relativity, the ADM one (more precisely, the equations arising 
from R a b = 0), and we choose exact lapse and exact co-shift (that is, the covariant shift). We use the same notation 
as in the body of the paper, except that now N = N(t,r) and (3 r — /3 r (t,r) are prescribed as arbitrary functions of 
spacetime. In spherically symmetry, the evolution equations for such a formulation are 
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2f3' r - 2NK r 



9t 



Mt +2 P I 9t_ 2NKt 



rg r 



K r 



(-2g^rf3 r K' rr g rr - Ag%rK rr {3' r g rr + Ag%rK rr (i r g' rr + 



K T 



2N 9rr9T r 9T + 4Ng 2 r g' T g T - Ng rr g' rr rg T g' T - 2N g rr g' rr g\, - Ng 2 r (g' T ) 2 r + 
2Ng rr K 2 „.g 2 T r - ANg 2 rr K rr g T rK T + 2N" g 2 rr g 2 T r - g' rr N' 'g rr g^ 
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(Bl) 
(B2) 

(B3) 

(B4) 
(B5) 

(B6) 
(B7) 



Note that in the previous equations, only gx appears with second derivatives. This means that we can rewrite 
this system as a first-order one by introducing just one new variable, z := g' T \ i.e. ii = Au' + B, where u — 
(g rr ,g T) K rr ,KT,zy. 

As an evolution equation for z we make the simplest possible choice: we just take the spatial derivative of the 
r.h.s. Eq. (B2) (i.e. we do not add the constraints to this new equation). Also, we replace everywhere g' T by z. The 
principal part is, then: 
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The determinant of the matrix that diagonalizes A is proportional to 

(-AK rr (3 r + g rr N'){rz + 2g T ) 2 
6AgU 2 l3 r r 2 



(B8) 

(B9) 
(B10) 

(BH) 
(B12) 



(B13) 



the system being strongly hyperbolic unless this determinant is zero. For the usual slicings of Schwarzschild (Painleve- 
Gullstrand, Kerr-Schild, time harmonic and full harmonic slicings) expression (B13) is different from zero, so the 
system is, indeed, strongly hyperbolic, and the same will hold for (perhaps slight) distortions of those spacetimes. 
Now, it is known that one has to give initial data and boundary conditions to the incoming characteristic modes 
in order to fix the solution in such a system. The point is that here there is always a negative eigenvalue (whose 
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corresponding eigenmode will thus propagate in the direction of increasing r, in particular, leaving the black hole): if 
[3 r 7^ 0, these eigenmode is either u± or 112, while the mode is U5 if (3 r = 0. 

It is important to point out that whether modes leave or not the black hole strongly depends on the particular 
formulation that one is dealing with. For example, using exactly the same formulation (ADM in spherical symmetry) 
with other choices of lapse and shift (such that they 'lock' the area) also gives strongly hyperbolic formulations when 
rewritten in first order form, but they do not have modes leaving the domain (see [|23|). Also, if one uses exact-lapse 
and exact-shift (as opposed to exact co-shift), the system is only weakly hyperbolic and, thus, the characteristic modes 
are not complete but, in any case, they do not leave the black hole (see, also, |p^| ). What we want to emphasize here 
is that superluminal modes can appear quite naturally, even in standard formulations. 
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